Two-dimensional nonlinear vector states in Bose-Einstein condensates 
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Two-dimensional (2D) vector matter waves in the form of soliton-vortex and vortex-vortex pairs 
are investigated for the case of attractive intracomponent interaction in two-component Bose- 
Einstein condensates. Both attractive and repulsive intercomponent interactions are considered. 
By means of a linear stability analysis we show that soliton-vortex pairs can be stable in some 
regions of parameters while vortex-vortex pairs turn out to be always unstable. The results are 
confirmed by direct numerical simulations of the 2D coupled Gross-Pitaevskii equations. 
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I. INTRODUCTION 

Multicomponent Bose-Einstein condensates (BECs) 
have been subject of growing interest in recent years 
as they open intriguing possibilities for a number of im- 
portant physical applications, including coherent storage 
and processing of optical fields [3, Q , quantum simula- 
tion Q, quantum interferometry etc. Experimentally, 
multicomponent BECs can be realized by simultaneous 
trapping of different species of atoms [j, 0] or atoms 
of the same isotope in different hyperfine states. Mag- 
netic trapping freezes spin dynamics 0, 0] , while in opti- 
cal dipolc traps all hyperfine states arc liberated (spinor 
BECs) Q . Theoretical models of multicomponent BECs 
in the mean-field approximation are formulated in the 
framework of coupled Gross-Pitaevskii (GP) equations 
and the order parameter of multicomponent BECs is 
described by a multicomponent vector. 

Like in the scalar condensate case, various types of 
nonlinear matter waves have been predicted in multi- 
component BECs. They include, in addition to ground- 
state solutions [13, EH, EH, structures which are pecu- 
liar to multicomponent BECs only, such as bound states 
of dark-bright [l3| and dark-dark [lj] , dark-gray, bright- 
gray, bright-antidark and dark-antida rk II 51 complexes of 
solitary waves, domain wall solitons jl6l. Il7l. ITc|. soliton 
molecules [l9j], symbiotic solitons [2(|. Two-dimensional 
(2D) and three-dimensional (3D) vector solitons and vor- 
tices have been considered in Refs. [2l|, [22|, [23[ for the 
case of repulsive condensates. Attractive intracomponent 
interaction have, on the other hand, received less at- 
tention and only one-dimensional vector structures have 
been studied so far [24], US]- Two-dimensional and 3D 
cases, however, demands special attention since the phe- 
nomenon of collapse is possible in attractive BECs. 

Interactions between the atoms in the same and dif- 
ferent states can be controlled (including changing the 
sign of the interactions) via a Feshbach resonance. The- 
oretical and experimental studies have shown that inter- 
component interaction plays a crucial role in dynamics of 
nonlinear structures in multicomponent BECs. Recently, 
two-component BECs with tunable inter-component in- 
teraction were realized experimentally [2(| [27|. Note 



that in nonlinear optics, where similar model equations 
(without the trapping potential) are used to describe the 
soliton- induced waveguides [28| , the nonlinear coefficients 
are always of the same sign. 

The aim of this paper is to study 2D nonlinear local- 
ized vector structures in the form of soliton-vortex and 
vortex-vortex pairs in a binary mixture of disc-shaped 
BECs with attractive intracomponent and attractive or 
repulsive intercomponent interactions. Then, by means 
of a linear stability analysis, we investigate the stability 
of these structures and show that pairs of soliton and 
single-charged vortex can be stable both for attractive 
and repulsive interactions between different components. 
Vortex-vortex pairs turn out to be always unstable. The 
results are confirmed by direct numerical simulations of 
the 2D coupled Gross-Pitaevskii equations. 

The paper is organized as follows. In Sec. [Til we for- 
mulate a model and present basic equations. The cases 
of attractive and repulsive intercomponent interactions 
are considered in Sees. Mil and IIVI respectively. The 
conclusions are made in Sec. IVl 



II. BASIC EQUATIONS 

We consider a binary mixture of BECs, consisting of 
two different spin states of the same isotope. We assume 
that the nonlinear interactions are weak relative to the 
confinement in the longitudinal (along z-axis) direction. 
In this case, the BEC is a "disk-shaped" one, and the GP 
equations take an effectively 2D form 
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where M is the mass of the atoms, 14 x t(r) 
y 2 )/2 is the harmonic external trapping potential with 
frequency lu±_ and V 2 = d 2 /dx 2 + d 2 /dy 2 is the 2D 
Laplacian. Atom-atom interactions are characterized 
by the coupling coefficients gij = Anh 2 aij/M, where 



2 



the s-wave scattering lengths for binary 
collisions between atoms in internal states \i) and \j). 
Note that gu — g 2 2 and 512 = 521- Introducing dimcn- 
sionless variables (x,y) — ► (x/l,y/l), t — > u±t, tyj — > 
^ 3 V?i^_ L /( 2|gii|), gip = — 512/l.gn I , -Bn = -gu/\gu\, 
where I = \JUj (MwjJ, one can rewrite equations |T]) and 
© as 
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In what follows we consider attractive interaction be- 
tween atoms of the same species and set B\\ = B22 = 1- 
We neglect the spin dynamics (assuming magnetic trap- 
ping) so that the interaction conserves the total number 
Nj (j — 1,2) of particles of each component 
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and energy 
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where 



Ej = y||V^| 2 -i|^| 4 + (x 2 + y 2 )|^| 2 }d 2 r. (7) 



III. ATTRACTIVE INTERCOMPONENT 
INTERACTION 

A. Stationary soliton-vortex pairs 

We look for stationary solutions of Eqs. © and © in 
the form 



(8) 



where rrij is the topological charge (vorticity) of the j-th 
component, fj,j is the chemical potential, r = \J x 2 + y 2 
and ip is the polar angle. Substituting Eq. ([5]) into Eqs. 
© and igj), we have 

/ir^i + A[ mi Vi ~ r 2 ^ + (l^il 2 + B 12 \tP2\ 2 )iPi = 0, (9) 
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FIG. 1: Normalized numbers of atoms N 1 and N2 of each 
component versus chemical potential /X2 at fixed ^ti for vector 
soliton-vortex pair (mi = 0, 7712 = 1) (attractive intercompo- 
nent interaction). 



1 1 1 ' 1 • r 



>--e--o--©--o--o--e--o--©--o--o--© 



STABLE 

a a— A ~~ y 



CM 



- -as* 9 * 



UNSTABLE 



-> — 1 — > — 1 — » — 1 

-5 -4 -3 -2 



~i • 1 ' 1 ' r~ 

-10 12 



FIG. 2: Existence domain for vector state mi = and m,2 = 1 
on the Qui, ^2) plane for attractive intercomponent interaction 
B12 = 0.5. Open circles correspond to numerically found 
existence boundaries. Dashed curves outline the variational 
predictions. At upper and lower boundaries of the existence 
domain the vector states degenerate into the pure scalar states 
with Ni — and ./V2 = respectively. The solid line with 
triangles indicates the stability boundary. 



^2 + A(™ 2 ) V2 - r 2 V2 + {Bu\ipi | 2 + |^2 1 2 )^2 = 0, (10) 

where Ar = d 2 /dr 2 + (l/r)(d/dr) — m 2 /r 2 . As was 
pointed out the inter-component interaction may be 
varied over wide range, however, the strength of the 
intercomponent interaction is weaker than the intra- 
componcnt counterpart in most experiments with two- 
component BECs. Thus -B12 can be considered as the 



free parameter from the range — 1 < -B12 < 1. We find 
that the qualitative behavior of vector state character- 
istics does not change when varying B12. To make it 
definite we further fixed the strength of intercomponent 
interaction at B\2 = ±0.5 for attractive and repulsive 
cases respectively. In this section we consider the case 
of attractive interatomic interaction B\2 > 0. First, 
we present a variational analysis. Stationary solutions of 
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FIG. 3: The numerically found solutions of Eqs. Q and (|10|l for attractive intercomponent interactions B12 > at fixed 
/ii = —2 for different values of \±2- 




FIG. 4: Typical dependence of maximum growth rates for 
L — 1 and L = 2 azimuthal modes as functions of [12 at 
fixed fj,i, (here fii — —2), B12 = 0.5. Note that the widest 
instability domain has L = 2 mode. 



Eqs. j9J) and JTOJ) in the form Eq. ([8]) realize the ex- 
ternum of the energy functional E under the fixed num- 
ber of particles N\ and iV 2 . We take trial functions ipj in 
the form 



$j( r ) = hj [ — 

Qj 



exp 



(11) 



which correspond to the localized state (mi, m 2 ) with 
vorticities mi and m 2 for |1) and |2) components respec- 
tively, a,j and hj are unknown parameters to be deter- 
mined by the variational procedure. The parameters hi 
and /12 can be excluded using the normalization condi- 
tions ((5|), which yield the relation Nj = nijlirh'ja'j. Thus, 
the only variational parameters are a\ and a^. Substitut- 
ing Eq. (fTTj) into Eq. (JT2J) , we get for the functional E 
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FIG. 5: Snapshots of unstable evolution for vector pair mi = 
0, m2 = 1, Hi = —2, [12 = 0.5, B12 — 0.5. The absolute 
values of \ipi\ and \ip2\ are shown in grayscale: the darker 
regions correspond to higher amplitudes. 
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2|mi| 2|m2| 



E U - ~2 B ^iiV 2(ffl?+a 2 )1+|mi|+|m2| 



(mi + m 2 )\ 



7rmi!m2! 



By solving the variational equations dE/daj = at fixed 
Nj one finds the parameters of approximate solutions 
with different mi and m 2 . We will focus on one par- 
ticular configuration with mi = and m 2 = 1 which 
corresponds to the pair soliton- vortex. The results of the 
variational analysis for this case and £?i 2 > are given 
in Fig. [T]and Fig. H]by dashed lines. These results were 
the starting point for numerical analysis. 

The equations ((9]) and (fT0|) were discretized on the 
equidistant radial grid and the resulting system was 
solved by the stabilized iterative procedure similar to 
that described in Ref. (2j|. The appropriate initial 
guesses were based on the variational approximation. 
The numerical results are shown in Fig. [1] and Fig. [2] 
by open circles, ft is seen that the variational results 
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exhibit a good agreement with numerical calculations. 

The stationary vector states form two-parameter fam- 
ily with parameters fii and /i2- In the Fig. [T]thc number 
of atoms for each component of the stationary vector 
state (0, 1) is represented as a function of the chemical 
potential /i2 at fixed value of /ii . The existence domain is 
bounded and its boundaries are determined by the con- 
dition ^2\m 1= q ^^2 5: M2|at 2= o- F° r each value of /ii, 
where the solution exists, a dependence similar to one 
presented in Fig. [T]can be found. This allows one to re- 
construct an existence domain of the vector pair (0, 1) on 
(/Xi, plane, which is shown in the Fig. [2] As is known, 
(see e.g. |3(|) for the two-dimensional scalar solitary 
structures in BEC with attraction, the chemical poten- 
tial is bounded from above fi < ^ max = 2 (to + 1), where 
to is the topological charge, and N — ► when fx — > /x max . 
One can see from Fig. [2] that the value of /x max is reduced 
in the presence of the second component if the intercom- 
ponent interaction is attractive (-B12 > 0). Both compo- 
nents vanish at the point (/xi, fj.2) = (2mi + 2, 2to2 + 2). 

Examples of soliton-vortex (0, 1) radial profiles are 
given in Fig. [3] It is interesting to note that when the 
amplitude of the vortex component is sufficiently high, 
a soliton profile develops a noticeable plateau. Such a 
deviation from gaussian-like shape leads to comparable 
divergence of numerical and variational dependencies in 
Fig. [Q Other vector states as (0,2), (-1,1), (-2,2) 
etc. were also found, but they all turn out to be always 
unstable (see below). 





FIG. 7: Numerically found solutions of Eqs. (|9} and (flu) for 
repulsive intercomponent interactions at fixed fj,\ = —2 from 
stable (right panel) and unstable (left panel) regions. 




FIG. 8: Same as in Fig. g] for B 12 = -0.5. Note that in 
contrast to the case of attractive intercomponent interaction 
the stability threshold is determined by L — 1 mode. 



states. We take the wave functions in the form 



^■(r,i) = {^-(r) +e j (r,t)}e- i ^ t+im ^, 



(12) 



FIG. 6: The same as in Fig. [2] for repulsive intercomponent 
interaction (B12 = —0.5). 



where the stationary solutions ipj{r) are perturbed by 
small perturbations £j(r, t), and linearize Eqs. ([3]) and 
(HI with respect to Ej. The basic idea of such a lin- 
ear stability analysis is to represent a perturbation as 
the superposition of the modes with different azimuthal 
symmetry. Since the perturbations are assumed to be 
small, stability of each linear mode can be studied inde- 
pendently. Presenting the perturbations in the form 

ej[r,t) = Uj (r)e luJt+lLv + v* j {r)e* , '* t ^ iI " p , 

we get the following linear eigenvalue problem 



B. Stability of stationary solutions 
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U = ujU, 



(13) 



The stability of the vector pairs can be investigated 
by the analysis of small perturbations of the stationary 



where U = (ui,Vi,U2,Va) is the vector eigenmodc and 
co is an (generally, complex) eigenvalue, aj = /3±2 = 



5 







1>1 



o • 1 



T=41.25 



T=41 .25 



2 I - 2 

I 

o P k-i 

-2 I -2 

T=|42.50 T=k2.50 



FIG. 9: Development of snake-type instability (L = 1) for 
vector pair mi = 0, rri2 = 1, /ii = —2, /i2 = 3.5. Intercom- 
ponent interaction is repulsive, B12 = —0.5. The absolute 
values of and \ip2\ are shown for different times. 
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integer i determines the number of the azimuthal mode. 
Nonzero imaginary parts in oj imply the instability of the 
state \ipi,ip2) with 7^ = max|Imoj| being the instability 
growth rate. 

Employing a finite difference approximation, we nu- 
merically solved the eigenvalue problem (TlU)) . Typical 
dependencies of the growth rate 7^ of the azimuthal per- 
turbation modes L = 1 and L — 2 on fX2 at fixed /xi are 
shown in Fig. 0]for the state (0, 1). Note that above some 
critical value of 1x2 ah growth rates vanish granting the 
stability of the vector pair (0, 1) against the azimuthal 
perturbations. For the case of attractive intercomponent 
interaction B12 > the stability boundary is determined 
by L = 2 mode. Similar dependencies of the growth rate 
7l on fi2 have been obtained for other values of fj,\. We 
have performed the numerical calculation of 7^ for val- 
ues of the azimuthal index up to L = 5. In all studied 
cases azimuthal stability is always defined by vanishing 
of the growth rate of L = 2 mode. The corresponding 
stability threshold is given in Fig. [2] by filled triangles. 
Note that in the degenerate scalar case ipi = ® for single- 
charge vortex ni2 = 1 the stability threshold 1x2 = 2.552 
coincides with the value obtained in Ref. [3fJ]. For the 
vector states (0, 2), ( — 1, 1), (—2, 2), the growth rates are 
found to be nonzero in the entire existence domain and 
these pairs appear to be always unstable. 

To verify the results of the linear analysis, we solved 
numerically the dynamical equations ([3]) and ^ initial- 
ized with our computed vector solutions. Numerical inte- 
gration was performed on the rectangular Cartesian grid 
with a resolution 512 2 by means of standard split-step 
fourier technique (for details sec e.g. 3 lj ) . In full agree- 
ment with the linear stability analysis, the states (0, 1) 



perturbed by the azimuthal perturbations with different 
L survive over huge times provided that the correspond- 
ing Hi and 1x2 belong to the stability region. On the 
other hand, Fig. [3] shows the temporal development of 
azimuthal L = 2 instability for the vector state (0, 1) 
with fii = —2 and /12 = 0.5 (i. e. in the instability re- 
gion). One can see the two humps which appear on the 
initially smooth ring-like intensity distribution. Further, 
the vortex profile is deformed, vortex and fundamental 
soliton both split into two filaments which then collapse. 
Note that unstable (0, 1) pair in BECs with repulsive in- 
tracomponent interaction [22j does not collapse and un- 
dergo a complex dynamics with trapping one component 
by another. 



IV. REPULSIVE INTERCOMPONENT 
INTERACTION 

In this section we present results for the case of repul- 
sive interactions B12 < between different components. 
The existence domain, stable and unstable branches on 
the (/j,i,/j,2) plane for the state (0,1) are shown in Fig. 
[5] It is seen that the repulsive intercomponent interac- 
tion leads to an increase of the maximum chemical po- 
tential fx max for each component compared to the case 
fx = 4 of pure scalar solution. The stability properties of 
the vector states were investigated by the linear stability 
analysis described in the preceding section. The states 
(0, 2), (—1, 1), (—2, 2) are always unstable as in the case 
of B 12 > 0. 

Figure [7] shows examples of the radial profiles of un- 
stable and stable (0, 1) states. In Fig. [5] we plot the 
growth rates 7^ as functions of 1x2 under fixed fi\ for the 
azimuthal perturbations with L = 1 and L = 2. The 
growth rates vanish if 1x2 exceeds a some critical value. 
In contrast to the attractive intercomponent interaction 
case, it is seen that the stability boundary is controlled 
by the elimination of the snake-type instability (i. e. az- 
imuthal perturbation with L = 1). Indeed, the repulsion 
between components naturally leads to spatial separation 
of condensate species. This relative motion destroys the 
vector state as seen from Fig. |H1 



V. CONCLUSIONS 

In conclusion, we have analyzed the stability of 2D 
vector matter waves in the form of soliton-vortex and 
vortex- vortex pairs in two-component Bosc-Einstcin con- 
densates with attractive interactions between atoms of 
the same species. Both attractive and repulsive intercom- 
ponent interactions are considered. We have performed 
a linear stability analysis and showed that, in both cases, 
only soliton-vortex pairs (0, 1) can be stable in some re- 
gions of parameters. Namely, under the fixed number of 
atoms in the soliton component, the number of atoms of 
the vortex component should be less than a some criti- 
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cal value. No stabilization regions have been found for firmed by direct numerical simulations of the 2D coupled 
vortex- vortex pairs and they turn out to be always un- Gross-Pitaevskii equations, 
stable. The results of the linear analysis have been con- 
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